Emergent microrobotic oscillators via asymmetry-induced order

Spontaneous oscillations on the order of several hertz are the drivers of many crucial processes in nature. From bacterial swimming to mammal gaits, converting static energy inputs into slowly oscillating power is key to the autonomy of organisms across scales. However, the fabrication of slow micrometre-scale oscillators remains a major roadblock towards fully-autonomous microrobots. Here, we study a low-frequency oscillator that emerges from a collective of active microparticles at the air-liquid interface of a hydrogen peroxide drop. Their interactions transduce ambient chemical energy into periodic mechanical motion and on-board electrical currents. Surprisingly, these oscillations persist at larger ensemble sizes only when a particle with modified reactivity is added to intentionally break permutation symmetry. We explain such emergent order through the discovery of a thermodynamic mechanism for asymmetry-induced order. The on-board power harvested from the stabilised oscillations enables the use of electronic components, which we demonstrate by cyclically and synchronously driving a microrobotic arm. This work highlights a new strategy for achieving low-frequency oscillations at the microscale, paving the way for future microrobotic autonomy.


Supplementary Notes 1 Mechanistic Model and Simulation of Particle Beating
As the microparticles beat at the curved liquid-air interface of a drop of aqueous H 2 O 2 solution, we start by solving the Laplace equation of capillarity [1,2]. We solved the following system of ordinary differential equations (ODEs) with MATLAB's ode45 Runge-Kutta solver (MathWorks, Inc., Natick, MA): ρ(0) = z(0) = θ(0) = V (0) = 0 (5) where β is the curvature at the apex s = 0 and V is the volume. γ c = g(ρ l − ρ a )/γ denotes the capillary constant, where γ is the interfacial tension, ρ l density of the liquid phase, and ρ a that of air. θ, r, z, and s are defined in Supplementary As discussed in the main text, the microparticles are driven outward by the collapse of a shared bubble and come together via a global and a local restorational force. As the SU-8 polymer is barely denser than the peroxide solution, the buoyancy from a small gas bubble underneath the disc is able to overcome the particle's weight and create a net force upward. As the microparticle is constrained to the liquid-air interface, it climbs the global drop profile defined by the above solution of the Laplace equation. One can formulate the energy as the product of the particle's vertical displacement and its weight after the subtraction of the Archimedes force. Thus, the lateral component of this global restorational force is given by: where V b (t) is the instantaneous bubble volume. Only the mass of the particle, m, is considered as that of the bubble is insignificant. The dimensionless factor Λ b is the volume fraction of the gas bubble lying below the undisturbed interface, as it is the displaced liquid in this region that gives rise to buoyancy. Note that we only included Λ b for the generality of Supplementary Equation (6). We use a Λ b of unity in the simulations hereafter in accordance with experimental observations. Given the drop profile, the force always points towards the apex. To quantify the "Cheerios effect", i.e. the inter-particle capillary attraction as a result of the local interfacial distortion, we adopt the Nicolson approximation [3] which assumes that (i) the horizontal force from capillary pressure is insignificant compared to that from buoyancy, and (ii) the small interfacial distortions may be superposed [4]. Prior results show that the Nicolson approximation is justified for small Bond numbers B = R 2 /L 2 c , or equivalently if the floating object's radius R ≪ L c = γ/ρ l g, the capillary length. Indeed, the L c of our experimental system is approximately 2.7mm, far exceeding the spatial scale of the beating physics. The surface height in the neighbourhood of a floating bubble follows: where l is the lateral distance from the bubble centre, K n the modified Bessel function of the second kind of order n, and Σ the buoyancy-corrected dimensionless weight defined by 2πγRBΣ = [m − ρ l Λ b V b (t)] g. Supplementary Equation (7) is a simplified asymptotic result true for l ≪ L c . The lateral capillary force experienced by a bubble of volume V ′ b at a distance l away is therefore: Needless to say, the capillary attraction force points towards the centre of the other particle. The K 1 (l/L c ) dependence is in agreement with the results derived from an energy approach [5]. Readers are directed to [6] for the treatment of scenarios with more than 2 particles. We next consider the hydrodynamic interactions. In the regime of low Reynolds number and low capillary number such as our system, the drag force for an object at the liquid-air interface is expressed as: where µ is the liquid's dynamic viscosity, R b the bubble radius, and v the instantaneous velocity. The drag coefficient Λ d is a scaling factor depending on the object's geometry, its depth of immersion, the contact angle, surface tension, and the densities of the object and the liquid [6]. As Λ d is difficult to estimate analytically, we assume it is a constant for simplicity and leave it as one of the two free parameters we estimate from experiments, a practice consistent with published models of microparticle motion along a curved interface [2]. An important additional consideration is the significantly increased drag when multiple particles approach one another, caused by the increased resistance to removing the liquid between them [7]. We note this inter-particle hydrodynamic interaction particularly because of the noticeable deceleration in our beating system when the edge-to-edge distance between particles were less than 2R p (see, for example, Fig. 1g between 68 and 69s). The approach velocity was virtually 0 right before contact, suggesting a drag significantly larger than that given by Supplementary Equation (9). Indeed, some previous studies predicted two floating microparticles to accelerate towards each other all the way until they collide if the Stokes' drag expression was not corrected for inter-particle interactions [4].
The most numerically convenient means of accounting for said interactions is to adopt the concept of hydrodynamic mobility [8], as with a number of previous works [6,9]. This correction factor as a function of the inter-particle spacing, l, is given by: where λ = l/ max[R b (t), R p ]. G, which is typically multiplied to the terminal velocity, approaches 1 for large separations λ − → ∞ and 0 for λ = 2 when the objects contact. Equivalently, we divided the drag expression in Supplementary Equation (9) with G in our numerical simulations.
The force expressions in Supplementary Equations (6), (8), (9), and (10) allow us to simulate the motion of each beating particle i with Newton's second law: We followed [2] in introducing a scaling factor for the effective mass (m eff = Λ m m) to account for the added mass of liquid experienced during particle acceleration. The two fitted parameters of the model, Λ m = 11.25 and Λ d = 0.35, were kept constant across simulations of different H 2 O 2 concentrations. The model outlined above was solved numerically again with MATLAB's ode45 Runge-Kutta solver. As all three forces are also dependent on the instantaneous bubble volume (equivalently, the radius), we zoom in to the catalytic surface and study the reaction kinetics as the final piece of the puzzle. The volume of the O 2 bubble as a function of time, V b(t) , is dictated by the rate of O 2 generation, which in turn is dependent on the free platinum patch surface area A Pt,free , as well as the peroxide concentration [H 2 O 2 ]. For a given experiment, we assume that the peroxide is in excess and [H 2 O 2 ] is a constant throughout, based on the absence of a shift in the beating frequency ( Supplementary Fig. 6). A well-studied catalytic reaction, the decomposition kinetics of H 2 O 2 on noble metal and oxide surfaces can be described by the classic Langmuir-Hinshelwood mechanism [10,11]: where k is a constant encompassing the reaction rate constant, the specific volume of O 2 , and the areal density of the surface sites. The kinetic equation represents that the rate is first order with respect to the concentration of bound surface sites, which are saturated at increasing peroxide concentration modulated by the binding constant K H . In the single particle scenario, A Pt,free decreases over time as the bubble underneath the particle starts to limit the accessible catalytic surface area. This leads to a reduced dV b /dt and therefore a self-limiting reaction. Inspection of the 2-particle beating videos, on the other hand, shows a near-linear increase of the bubble volume up until the moment of merger. This observation suggests that the bubbles in the beating system do not grow beyond the critical V b which marks the onset of catalytic surface blockage, and that A Pt,free ≈ A Pt .
Vertical coordinate relative to the drop's apex [m] 2 Thermodynamics of Particle Beating Asymmetry-induced order is a process by which explicit symmetry-breaking (spontaneous or otherwise) leads to the emergence of ordered states in a system [12][13][14][15]. Hence, asymmetry-induced order requires both a symmetry whose breaking can be observed, and a clear notion of "degree of order." Which symmetry to break is inherently a system-dependent question, and as such there are no general means of choosing between symmetry groups to achieve a desired outcome. However, the so-called degree of order of a system is a challenging property to formally specify in general. For one, what is meant by order is often ill-defined or underspecified. Secondly, even when provided with a means to metricize order, such metrics are often analytically and computationally intractable because they require global knowledge of system states-as is the case for calculating entropy. This is further complicated by the fact that, far from equilibrium, entropy is not sufficient to establish the robustness, stability, or persistence of system configurations (all of which are attributes often ascribed to "orderly" states) [16]. To this end, physicists have made use of order parameters to establish more narrowly-construed notions of order on a case-by-case basis for particular systems [17,18]. Recent work in nonequilibrium thermodynamics has made strides towards describing the emergence of order more generally in broader classes of complex systems. Rattling theory is a novel thermodynamic theory describing the emergence of order and self-organization in "messy" nonequilibrium dynamical systems [19,20]. The success of rattling theory depends crucially on the definition of the class of systems it considers to be messy. The rattling ansatz sees the behavior of complex systems as stochastic diffusion processes taking place in high-dimensional configuration spaces in the presence of energy influxes driving them out of equilibrium. Any system whose behavior can be described by such configuration-space diffusion falls under the class of messy systems described by rattling theory. Modelling the behavior of systems as diffusion processes is what enables an analytical determination of nonequilibrium steady-state density, and, as a consequence, an understanding of self-organization. Empirically, this approach has been shown to predict the long-term behavior of a wide variety of systems, from canonical chaotic systems [21] to swarms of robots [19], and is expected to apply across diverse active matter systems as well [22,23].
At the heart of the theory lies a precise, local, and computable measure of order, rattling, from which the theory derives its name. Rattling measures the way in which system configurations respond to external force fluctuations: Rapid, uncorrelated configurational changes produce high rattling values, and slow, correlated changes produce low rattling values. When a system's response to local force fluctuations is random (i.e., has Gaussian statistics), rattling is exactly the entropy of its configurational velocities. As a quantity, the rattling R(q) of a system at configuration q is where ⟨·, ·⟩ is the covariance tensor of the system's configurational velocities (i.e., two-point correlation function) averaged over an ensemble of dynamical trajectories initialized at q. Using this definition, we can state the central prediction of rattling theory, known as the "low-rattling selection principle." The principle expresses a relationship between the magnitude of system-level fluctuations measured at a particular configuration (i.e., rattling R(q)), and said configuration's prevalence in the system's nonequilibrium steady-state density. In particular, this relationship is of Boltzmann-like form: where γ is a constant of order 1. This relationship shows that configurations with remarkably low entropy dynamical responses (i.e., low rattling) are exponentially preferred in the system's steady-state. Thus, rattling and its associated selection principle are able to sufficiently establish the robustness, stability, and persistence of the configurations of complex systems. In summary, rattling captures the ways in which correlations among disorderly degrees of freedom give rise to system-level fluctuations of different magnitudes. Then, the low-rattling selection principle states that such system-level fluctuations bias the nonequilibrium steady-state of a complex system towards configurations in which the system experiences remarkably low magnitude fluctuations. Furthermore, this spontaneous selection of low rattling configurations necessarily requires that strong correlations between degrees of freedom arise, and thus for orderly behaviors to emerge. Interested readers looking for a complete treatment of this material, as well as theoretical derivations and experimental validation, are referred to [19].
Equipped with a precise way to quantify order in a broad class of complex systems, we may now develop a system-specific understanding of the ways in which symmetry-breaking affects the rattling of our system of beating particles in hopes of finding strategies to stabilize periodic system beating for N > 2.
In order to elucidate the role that symmetry-breaking may play in the self-organized states of our system of active microparticles, we must now consider specific system symmetries and their relationship to the magnitude of system-level fluctuations. While our system is not invariant to the action of any obvious continuous symmetry groups, it is permutationsymmetric [24]. This is to say that our collection of microparticles are all dynamically identical (up to fabrication tolerances). Hence, one promising avenue to investigate is the different ways in which permutation-symmetry breaking may lead to order in our system. Based on results from our mechanistic modelling of particle beating, we know that there are two ways in which the dynamics of individual microparticles can be made distinct from one another. First, we know that changing the volumetric shape of particles will lead to different local hydrodrynamic drag properties. Second, we know that changing the buoyancy of particles also produces local changes to individual microparticle dynamics through its effect on capillary forces. However, changing the shape of our microparticles requires major changes to their fabrication, as well as nontrivial modifications to the mechanistic model. In contrast, we can easily modify a particle's buoyancy by modulating the volume of the bubble forming underneath the particle, which we can in turn control through the size of their Pt patch.
To explore the role of permutation-symmetry breaking on our system, we constructed a simple model that we can work with analytically from the perspective of rattling theory. In line with the rattling ansatz, our model considers the configurational dynamics of collectives of beating particles as a diffusion process. We incorporate the effect of heterogenous particle buoyancies through the inclusion of a parameter modulating the size of bubbles in analogy to the role of the Pt patch. Our beating particles are perfectly suited for this sort of analysis, even more so than others (e.g., robot swarms). In part, this is due to the physics of fluid dynamics at low-Reynolds numbers (∼0.25 Re for our system) [25]. In this regime, inertia ceases to influence the behavior of systems, leaving viscous forces and stochastic thermal fluctuations to affect their dynamics substantially-thereby making a diffusive approximation natural.
Our model elucidates the role of design parameters on the structure of the system-level fluctuations on the basis of two primary assumptions. First, we assume that the behavior of each individual particle i is monotonically modulated by some real-valued design parameter U i from a set U = {U 1 , · · · , U N } for a system of N particles. These design parameters correspond by analogy to the Pt patch size. Second, we assumed that particle i's bubble burst only affects the other members of the collective and not itself, which broadly matches experimental observations. We can think of the U i parameters as implicitly determining the strength of the impulse imparted by particle i's bubble burst onto its neighbors. In particular, we model the effect of this parameter and the bubble burst strength a i according to the following Boltzmann-like monotonic relationship, where Z is a normalization factor given by Z = N i=1 e −Ui . In other words, the a i parameters can be thought of in analogy to the size (and strength) of bubbles that a given particle can support. Hence, we can motivate this modeling choice by envisioning the gas in bubbles distributing itself according to an energy landscape specified by our U i parameters, and thusly influencing the bubble popping strength a i . The normalization factor Z arises from the fact that we are not interested in the absolute magnitude of the bubble bursts but rather the effect of their relative magnitudes on the collective behavior.
In the main text, we made use of an observable termed the "breathing radius" for purposes of analysis. The breathing radius is the mean Euclidean distance of the particles to the centroid of the collective. Similarly, here we will only consider the statistical properties of the dynamics of a breathing-radius-like observable,r(t), under a simple diffusive model. As in the main text,r(t) is an averaged quantity over particles:r(t) = N i=1 r i (t). By assumption, a bubble burst at particle i leaves particle i stationary, but a burst from some neighbor j exerts an impulse of random direction onto particle i. In this case, the dynamics of r i (t) evolve according toṙ where ξ j is normally-distributed delta-correlated multiplicative noise in the Itô convention. Note that this construction results in an anisotropic diffusion tensor without spatial dependence, as we are not modelling the geometry of interparticle interactions but rather their statistical fluctuations. From this specification of the system's diffusive dynamics, we can apply rattling theory to understand the effect of our design parameters U i on the self-organized collective behavior of the system. Given this formulation of the system dynamics, we proceed by calculating the effect of parameter changes on the magnitude of system-level fluctuations. Letting r(t) = [r 1 (t), · · · , r N (t)] T , the correlation structure of the system is where δ kl is the Kronecker delta. We note that the correlation structure of the system has no dependence on time (i.e., it has infinite temporal correlation) and no dependence on configuration r(t), leaving the design parameters U i as the only variables with an effect on the system behavior. Finally, in order to express the system's rattling in terms of its design parameters we require an analytical expression for the determinant of its covariance tensor, which is challenging in general. Fortunately, for this particular correlation structure there exists such a closed-form expression, which enables us determine the system's rattling as a function of its parameters: Equipped with an understanding of how the system's parameters affect its rattling (and thus its degree of order), we can now use the model as a tool to guide our experimental design. While there are infinitely many parameter combinations for a Supplementary Figure 2: Rattling as a function of patch size in diffusive model. Here, we study the effect of a given particle's U parameter (in analogy to Pt patch size) on the rattling of collectives of varying sizes. Note that we subtract the constant offset in rattling due to system size so that R = 0 at U = 0 for all N . We find that any variability in the size of the particle's patch produces a drop in rattling, leading to asymmetry-induced order. When a particle becomes inert as U increases, it stops contributing to system-level fluctuations, leading to a modest drop in rattling independent of N . However, as U decreases the modified particle's bubble bursts dominate and effectively become the sole source of variance in the system's configurational degrees of freedom. Such coordination among degrees of freedom leads to a sharp drop in rattling dependent on N .
given collection of N particles, one of the simplest design alterations to study is the effect of a single particle differing from the rest-for reasons that we will see shortly, we term this particle a designated leader. In this setting, one particle will have its parameter be U DL while the rest of the N − 1 particles will have it beŪ (which we take to be a constant fixed a priori). (18), we have the following expression

Rearranging the expression in Supplementary Equation
which allows us to make predictions about the behavior of a collection of N beating particles with a single designated leader. However, much in the same way that entropy can trivially depend on system size (e.g., number of microstates), our expression for rattling in Supplementary Equation (2) does as well. Thus, to focus on the dependence of R(U DL , N ) on U DL , we subtract the constant bias that system size contributes to the value of rattling. To do this, we calculate R(U DL , N ) − R(Ū , N ) for a choice ofŪ that we fix across all system sizes, where we note that R(Ū , N ) is merely a constant that offsets the value of rattling to be zero when U DL =Ū . Since R(Ū , N ) is exclusively a function of the number of particles for a givenŪ , subtracting it from R(U DL , N ) precisely removes the constant contribution of system size to the overall magnitude of rattling. As detailed in [19], constant offsets to the rattling values of a system do not affect its behavior. Only changes to the rattling landscape-that is, changes to the relative rattling values between configurations (or parameters)-have an effect on system behavior. This implies that comparing the absolute rattling values across systems is of limited use, which motivates our approach (as in Fig. 2d and Supplementary Fig. 2).
In Supplementary Fig. 2 we show the results of varying the parameters of the designated leader for collectives of various sizes, while fixingŪ = 0 and subtracting the bias in rattling due to system size. Crucially, we observe that any deviation from the parameter values of the rest of members of the collective (i.e., away from U DL = 0) results in a reduction in rattling. Thus, our model predicts that any amount of heterogeneity will lead to increasingly ordered system states. Such asymmetry-induced order has been studied in networked systems of oscillators [12][13][14][15], but its emergence as a low-rattling phenomenon is a novel finding.
Through this mechanism, order arises in one of two distinct ways. First, as U DL increases, the designated leader particle becomes effectively inert. This is to say that the strength of its bubble bursts a DL asymptotically approach zero, as though it were a patchless particle. As a result, the leader particle acts as dead weight and does not contribute to system-level fluctuations, leading to a modest decrease in rattling-independent of the total number of particles-that matches experimental observations. Second, as U DL decreases, the designated leader particle's bubble bursts become stronger and its contribution to the magnitude of system-level fluctuations dominates over those of other particles. In turn, this effectively leads to a concentration of all variability and randomness in the system into a single of its many degrees of freedom, thereby leading to significant correlations in the behavior of all particles and a resulting drop in rattling. Note that as more particles are added more degrees of freedom become correlated, leading to sharper drops in rattling as a function of N . Hence, on the basis Supplementary Figure 3: Effect of designated leader on self-organization. By introducing a designated leader, the entropy of the bubble burst forcing patterns decreases (since they become periodic), which has an effect on the self-organization of the system. On the left panel, we simulate the dynamics in Supplementary Equation (16) and calculate their rattling and steady-state densities numerically. On the right panel, we consider experimental data from an 8 particle collective in both standard (∆U DL = 0%) and designated leader configurations (∆U DL = 40%), which we then process using the same procedure as for the left panel. While the absolute magnitudes of parameter values for the simulation are arbitrary, the ∆U DL values are determined from the actual Pt patch sizes used on the experimental systems. For both the simulated and the experimental data, the results are consistent with rattling theory (in particular, Supplementary Equation (19) with γ = 1). This procedure was then repeated for experimental collectives of other sizes with the same results. Hence, bubble burst patterns of varying entropy (which depend on system design parameters) provide an explanation for the emergence of system order that is consistent with our results. of these results and other studies of asymmetry-induced order we chose to study the influence of designated leaders on the collective behavior experimentally by producing leader particles with larger Pt patches.
Another consequence of applying the rattling ansatz to the behavior of our collective of beating particles is that it allows one to reinterpret the relationship between system elements. In the theory, configurations with exceptionally orderly responses to external driving forces are selected for the nonequilibrium steady-state of complex systems. Similarly, we can think of the relationship between the particle configurations (i.e., their relative positions and orientations) and the sequence of forces the system experiences due to bubble bursts according to this dichotomy. Working from this perspective, rattling theory then suggests that the entropy of the sequence of bubble bursts can prevent the system from finding orderly configurations by increasing their rattling (see [19], and in particular figure 4 within). More precisely, we have where S(q) is the entropy of the driving forces affecting the system at configuration q. Importantly, this expression shows that the effect of drive entropy is to simply offset a configuration's rattling.
In the main text, we observed that our system design parameters (i.e., the particle Pt patch sizes) do have a profound effect on system behavior and also on the bubble burst sequence-changing its behavior from seemingly random to almost perfectly periodic (see Fig. 2). If we were to accept the hypothesis presented by Supplementary Equation (19), then this difference in behavior should be explained by the difference in the entropy generated by the bubble burst patterns at different system parameters. Moreover, if this is the case, then the results from analyzing standard and designated leader systems should lie on the log-linear correlation of Supplementary Equation (19), with a slope of γ (which nominally is of order 1). Indeed, this is precisely what we observe in Supplementary Fig. 3 for simulations of the model dynamics and for our experimental data samples.
While throughout this section we have motivated the analytical model in specific reference to our experimental beating particles, our model and results generalize beyond our system. At its core, our generic model describes the structure of statistical fluctuations in a collection of strongly interacting degrees of freedom (i.e., under strong mixing conditions). This is to say that the details of how the magnitudes of said fluctuations are parametrized by system properties are inessential to the results. Particularly, the equation for rattling in Supplementary Equation (18) can be expressed in terms of a i directly as log (N − 1) N i=1 a i , which allows one to freely model the way in which individual degrees of freedom contribute to the overall system-level fluctuations. Thus, rattling theory, as well as the mechanism we have outlined for asymmetry-induced order, present a general framework from which to understand the effect of system design parameters on the self-organized behaviors of the system-providing a novel approach to micro-system design based on thermodynamic principles relevant to the scales of interest.

Note on Microelectronic Low-Frequency Oscillators
In this section, we elaborate on the design and fabrication challenges of microelectronic oscillators with a frequency on the order of a hertz, which we briefly alluded to in the Introduction part of the main text. Given the relatively large footprints of integrated capacitors and inductors available, RC-and LC-based oscillators are hardly compatible with the limited space on micrometre-sized machines [26]. For example, the frequencies of RC oscillators, such as a bi-inverter or a Schmidt Trigger oscillator, are on the order of the reciprocal of their respective RC constants, i.e. f RC ∼ O(1/RC). Taking the capacitance to be a generous 40pF for an area of 100µm×100µm [27], one would require a massive resistor of 25GΩ to achieve an RC time constant of 1s. Assuming a resistivity of 100kΩ/µm 2 of polysilicon, this resistor alone would occupy 2.5×10 5 µm 2 . Alternatively, one may opt to use a frequency divider to bring the kHz-order frequency of a typical microelectronic relaxation oscillator down to 1Hz. Suppose the starting frequency is 17kHz [28], a cascade of 15 T flip-flops is needed, each of which is constructed from at least 20 transistors [29]. Should 300 transistors be fabricated onto a 100µm×100µm microchip, the appropriate transistor node would be 500nm. While well within the realm of possibility, such technology typically still requires the involvement of a commercial foundry outside of academic institutions. Similarly, thyristor-based oscillators of frequencies from 20Hz and up have been foundry-fabricated with a feature size of 180nm [26]. The integrated circuit design expertise and capital investment required are the reasons for a high barrier-to-entry. Note that the area reserved for onboard energy harvesting and storage units, as well as for miscellaneous electronics, may further constrain the real estate available to the microelectronic oscillator.

Note on the Fuel Cell's Open-Circuit Voltage
Supplementary Fig. 13 shows that the open-circuit voltages of the Pt-Au and Pt-Ru fuel cell devices, V OC , exhibits a very weak dependence on the peroxide concentration [H 2 O 2 ], unlike the trend of the short-circuit current densities ( Fig. 4c and Supplementary Fig. 14). Here we provide a simple explanation based on electrochemical kinetics. We consider the following two pairs of forward and reverse reactions taking place on a single electrode: where ϕ • eq denotes the standard equilibrium potentials. The Butler-Volmer equation suggests that only one half-reaction from R 1 and R 2 each is dominant at the mixed potential ϕ mix , defined as the potential where the total current equals 0 [30]. If we consider the oxidative half-reaction of R 1 and the reductive half-reaction of R 2 (choosing the other two half-reactions does not alter the conclusion), the full Butler-Volmer kinetic expression is given by [31]: where ϕ is the applied potential on the absolute scale, i 1 (ϕ) and i 2 (ϕ) the respective current densities, n the number of electrons transferred, F the Faraday constant, k the rate constants, ν the reaction orders, α the the transfer coefficients,R the universal gas constant, T the absolute temperature. We can obtain the mixed potential ϕ mix by solving: which is equivalent to −i 1 (ϕ mix )/i 2 (ϕ mix ) = 1. While the exact form of the solution is of little relevance to us, the division of the right-hand side of Supplementary Equation (20)  5 Note on the Energy Expenditure

Energy Conversions of the Mechanical Oscillation
Within each period of the emergent mechanical oscillation, chemical energy stored in the H 2 O 2 fuel is converted into the particles' kinetic energy upon the collapse of the O 2 bubble. The kinetic energy imparted to two outgoing particles simply take the form of E k = mv 2 , where m is the mass of each particle and v the maximal velocity right following the bubble collapse. With m = 2.34µg for a 500µm-diameter particle and v = 3.2 × 10 4 µm/s measured from experiments, E out is estimated to be 2.40 × 10 −12 J per cycle.
The chemical energy consumed per cycle may be computed as: where V b,th denotes the bubble volume at threshold, estimated to be 9.81×10 −2 µL in a 2-particle homogeneous system in 1mL of 10% H 2 O 2 . We assume an ambient pressure P of 1atm and temperature T of 25 • C, as the excess Laplace pressure within the bubble before collapse is a negligible 5.0 × 10 −3 atm. ∆H, the enthalpy change of the decomposition reaction, is 98.24kJ/mol at given conditions, equivalent to an energy density of 2.89kJ/g H 2 O 2 or 0.29kJ/g 10wt% H 2 O 2 solution [33]. E chem per cycle is computed to be 7.88 × 10 −4 J. The portion of the chemical energy converted to the work of expansion is: where R b,th is the threshold radius assuming a spherical bubble. The latter term of 7.41 × 10 −8 J is the surface energy E surf , i.e. the work against the Laplace pressure during bubble growth. To summarize, therefore, 1.26% of the original chemical energy contributes to a W P V of 1.00 × 10 −5 J. 0.74% of the work of expansion is stored as the surface energy. Finally, the kinetic energy gained by the particles account for 0.032‰ of surface energy stored in the bubble.

Energy Conversions of the Microgenerators
As the microgenerator converts the chemical energy from H 2 O 2 decomposition to electrical work, it is of interest to calculate the proportion of total H 2 O 2 molecules consumed which contributed to the electrical current [34,35]. Given that each electrochemically redoxed H 2 O 2 molecule transfers an electron, ON-state currents of 180.66nA (in the absence of an electrical load) and 15.27nA (with a load, i.e. the actuator) are respectively attributed to 1.87×10 −12 and 1.58×10 −13 moles of H 2 O 2 per second. These correspond to 0.76‰ and 0.063‰ of the total peroxide consumption rate (2P/RT · dV b /dt = 2.45 × 10 −9 mol/s), respectively. The former is in agreement with prior literature [34], which estimated an electrochemical contribution of 0.5‰. Since more than 99.9% of the consumed H 2 O 2 decompose via the same non-electrochemical pathway as in the beating particles with no fuel cells aboard, generation of the electrical current has a negligible impact on the mechanical oscillation if all other conditions are kept the same. Along the same lines, additional fuel cell particles are not expected to diminish the electrical signals observed.

Supplementary Figures
Supplementary Figure 4: Beating particle fabrication steps. a, An array of SU-8 polymeric microdiscs were defined and patterned on a Si wafer with standard photolithography, followed by electron-beam physical vapor deposition of Pt on top. b, The particles were subsequently lifted off in heated KOH solution which etched into the Si substrate. The KOH was displaced by water in which the lifted off microparticles were stored. c, Alternatively, a film of PMMA polymer was spun over the microparticle array. Together they would delaminate from the substrate in heated KOH solution. The PMMA was then removed with an acetone rinse. The lifted off particles were transferred to water for storage.
Supplementary Figure 5: Detailed comparison between experimental and simulated beating behaviours of two particles. a, Mechanistic model simulations and experimental snapshots taken at representative stages of a beating cycle. b, The simulation and experiments are in excellent agreement, evident from the matching curves of the breathing radius, previously also shown in Fig. 1g. We note that the mechanistic model captures fine details of the self-oscillation, such as the subtle step change in (b) at approximately 69s. The step increase was a result of the merged bubble pushing the particles outwards slightly, reflected by both the experiment and simulation in (a). c, This panel shows the excellent agreement between the experimental bubble radii and those predicted by the mechanistic model. The former were measured manually from the raw video data.
Supplementary Figure 6: Beating frequencies over time from moving window recurrence analyses. The same histograms as in Fig. 1i of the main text were generated, but here only for breathing radius data within a moving window of 150 frames (5s). Frequencies calculated from the most probable recurrence time of each window were plotted as a function of time. The beating frequencies in all experiments are constant throughout, demonstrating robust periodicity. Furthermore, curves from experimental replicates overlap. The frequencies from moving window analyses agree with those shown in Fig. 1j for all the H 2 O 2 concentrations. These concentrations respectively correspond to 6-, 3-, 2-, and 1-fold volumetric dilution of a 30wt% H 2 O 2 solution.
Supplementary Figure 9: Robustness of the emergent oscillation to perturbations. In these two experiments, we intentionally disturbed a system of two identical particles by (i) deforming the liquid-air interface with a pipette [36], (ii) stirring the H 2 O 2 drop, and (iii) shaking the drop back and forth. It is evident in the breathing radius trajectories that the collective oscillation resumes promptly following the perturbations (shaded region) with its amplitude and periodicity unchanged, thus demonstrating robustness. Data discontinuities during the perturbations are a result of blurry frames or particles temporarily exiting the camera field-of-view. The inset micrograph shows the particles approaching the pipette due to the deformed interface. Scale bar, 1mm. Supplementary Figure 14: Short-circuit current density as a function of H 2 O 2 concentration for a Pt-Au device.
(cf. Fig. 4c of the main text). Error bar, standard deviation.
Supplementary Figure 15: Bimorph actuator experimental setup. The microactuator in PBS solution is connected via external wiring to the beating system in an H 2 O 2 drop. The mechanical self-oscillation is translated to an oscillatory electrical current as illustrated in Fig. 4a, which powers cyclic motion of the actuator (Fig. 4e).
Supplementary Figure 16: Oscillatory mechanical beating drives on-board oscillatory current. (See also Fig. 4e of the main text). As a standard 500-µm particle beats with a Pt-Ru fuel cell device (Fig. 4b, also Methods), the bubbles collapse at regular intervals as indicated by the spikes in the breathing radius trajectory (r(t), top). Removal of the bubbles restores the electrochemical reactivity of the fuel cell electrodes, and therefore the current (bottom) peaks precisely as r(t) does. The current measured in this experiment is an order of magnitude higher than that in Fig. 4e since the system characterized here was not connected to an actuator.